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Abstract 



We analyze the dynamics of a solid-state laser driven by an injected sinu- 
soidal field. For this type of laser, the cavity round-trip time is much shorter 
than its fluorescence time, yielding a dimensionless ratio of time scales o" <C 1. 
Analytical criteria are derived for the existence, stability, and bifurcations of 
phase-locked states. We find three distinct unlocking mechanisms. First, if the 
dimensionless detuning A and injection strength k are small in the sense that 
k = 0{A) <C o"^/^, unlocking occurs by a saddle- node infinite-period bifurca- 
tion. This is the classic unlocking mechanism governed by the Adler equation: 
after unlocking occurs, the phases of the drive and the laser drift apart mono- 
tonically. The second mechanism occurs if the detuning and the drive strength 
are large: k = 0(A) » o"^/^. In this regime, unlocking is caused instead by a 
supercritical Hopf bifurcation, leading first to phase trapping and only then 
to phase drift as the drive is decreased. The third and most interesting mech- 
anism occurs in the distinguished intermediate regime k,A = 0{a^/^). Here 
the system exhibits complicated, but nonchaotic, behavior. Furthermore, as 
the drive decreases below the unlocking threshold, numerical simulations pre- 
dict a novel self-similar sequence of bifurcations whose details are not yet 
understood. 
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I. INTRODUCTION 



The Adler equation 

^ = A-A;sin$ (1) 
dt ^ ' 

provides the simplest model of phase locking between a nonlinear oscillator and an external 
periodic drive. Here is the phase difference between the oscillator and the drive, A is the 
frequency detuning, and k is the coupling strength. This equation first arose in connection 
with the phase locking of microwave oscillators and has since found application in many 
other settings, including the depinning of charge-density waves 0, the entrainment of bio- 
logical oscillators P,^, and the onset of resistance in superconducting Josephson junctions 

ii- 

A system governed by the Adler equation can display only two types of long-term behav- 
ior 0. If |A/A;| < 1, all solutions tend to a phase-locked state, where the response oscillator 
maintains a constant phase difference relative to the driver. On the other hand, if |A/fc| > 1, 
all solutions exhibit phase drift, where the phase difference grows monotonically, with one 
oscillator periodically overtaking the other. 

The main limitation of the Adler equation is that it treats the response oscillator as 
a system with only one degree of freedom, namely its phase. Possible variations in its 
amplitude (and any other degrees of freedom) are ignored. This approximation is reasonable 
in the limit of weak driving; in that case, the amplitude of the response oscillator typically 
equilibrates much more rapidly than its phase, and can therefore be treated as a constant in 
the subsequent analysis. But if the driving is not weak (in some appropriate dimensionless 
sense), the dynamics can become complicated. In this paper we revisit a classic problem - 
the mathematical analysis of a solid-state laser with external injection [[^-[TT1| - and explore 



it in regimes where amplitude effects become important and the Adler approximation breaks 
down. 

Our work was inspired by recent theoretical and experimental studies of amplitude effects 
in two mutually coupled solid-state Nd:YAG lasers ||T^JT3[| . In those studies, the lasers 



were equally coupled and identical, except for a slight relative detuning of their frequencies 
from some common cavity mode. For coupling strengths well above or below the locking 
threshold, the lasers were found to exhibit the simple behavior expected from the Adler 
approximation. However, as the coupling approached the locking threshold from below, 
the lasers showed a series of amplitude instabilities, culminating in a period-doubling route 
to chaos. These novel instabilities could not be explained by the Adler approximation. 
Instead the authors proposed the following mechanism. Below the locking threshold, the 
lasers exhibit phase drift. If the time required for one full cycle of phase slip happens to 
be an integer multiple of the lasers' relaxation period, the resulting subharmonic resonance 
might account for the observed instabilities. For the highly symmetrical case where the two 
lasers are assumed to have identical intensities and gains, this argument was proven to be 
correct by reducing the governing equations to those for a single, periodically modulated 
laser, where the subharmonic resonance mechanism was already known to occur |jl^,|15|. 



We wondered whether similar amplitude instabilities and chaos would occur in two cou- 
pled Nd:YAG lasers with unidirectional coupling (or equivalently, in a single Nd:YAG laser 
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with external injection). On the one hand, the quahtative argument about subharmonic 
resonances should still work. On the other hand, the equally-coupled case enjoys special 
symmetries that are not present in the unidirectional case. Given the crucial role of the 



symmetry in the earlier analysis of Erneux et al. ||T2|,|13[, it seemed possible that some new 
effects might occur if the symmetry were broken. 

As we will show below, the system with one-way coupling can indeed display some 
fascinating behavior near the unlocking threshold, but it differs from that seen in the equally- 
coupled case. In particular, we do not see a period-doubling route to chaos, nor any evidence 
of chaos at all. Instead, in a certain distinguished regime of parameters, we find a self-similar 
cascade of periodic windows and bifurcations. To the best of our knowledge, this bifurcation 
scenario is novel. 

It will be interesting to see whether this cascade can be detected experimentally for a 
laser in the appropriate parameter regime, as specified by our theory. It would also be 
gratifying to have a better theoretical understanding of the cascade itself. 

This paper is organized as follows. The governing equations are given in Sec. In Sec. 
pn| , we reduce the number of parameters by nondimensionalizing the equations and exploiting 
certain symmetries. By choosing a frame that co-rotates with the phase of the driver, we 
reduce the system to three coupled autonomous ordinary differential equations: one for 
the phase difference between the laser and the drive, and one each for the dimensionless 
gain and amplitude of the response laser. Fixed points of this reduced system correspond 
to injection-locked states of the original system. Sec. |V| dispenses with the limiting cases 
of zero coupling or zero detuning where the dynamics can be analyzed straightforwardly. 

The analysis begins in earnest in Sec. |V|, where we derive criteria for the existence and 
stabihty of injection-locked states, and compare our criteria to those obtained in the usual 
Adler approximation. In Sees. |V| and |V11| , we start to investigate what happens when 
locking is lost. We show perturbatively that for a broad range of parameters, the phase 
difference oscillates periodically, but the precise nature of those oscillations depends on 
the relative sizes of the dimensionless coupling, detuning, and stiffness of the system. 

For a distinguished limit of parameters, described in Sec. VTT|, the reduced system 
has complicated dynamics and undergoes the self-similar cascade of periodic windows and 



bifurcations mentioned above. In Sec. IX, we consider the system in the singular limit of 



zero stiffness. Again, the cascade persists. Based on the distinctive helical structure of the 
periodic orbits, we propose a mathematical mechanism underlying the cascade. It apparently 
stems from a co dimension-two bifurcation in which a supercritical Hopf bifurcation combines 
with a saddle-node infinite-period global bifurcation. We have no proof of this mechanism, 
but show that it correctly predicts the scaling laws found numerically in the bifurcation 
diagram. We conclude with a discussion of open questions. 



II. FORMULATION OF THE MODEL 

For solid-state lasers, as well as other Class B lasers with negligible linewidth enhance- 
ment factors, such as CO2 and ruby (NMR) lasers p|, the polarization relaxes rapidly 
compared to the electric field and the gain, and can therefore be adiabatically eliminated. 



Following Ref. [IG] with straightforward modifications, we can write the following equations 
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for two lasers coupled through transverse overlap of their electric fields, assuming single- 
mode operation and neglecting spatial variations within the lasers: 



dT* 

rf,\pi-Gi-Gi\E 



dT* 
dE. 
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dT* 

^^2 _ ^-i,^ ^ \ -p \2 

dT* 



T, 



■f2\P2-G2-G2\E2\'). (2) 



Here T* is time, and for j = 1,2, Ej is the complex electric field, Gj is the gain, Tcj is 
the cavity round-trip time, Tfj is the fluorescence time of the lasing ions, pj is the pump 
coefficient, aj is the cavity loss coefficient, Uj is the detuning of the laser from some common 
cavity mode, and is a complex coupling coefficient with KEj representing the degree of 
overlap of the two lasers, with possible attenuation and dispersion taken into account. The 
noise term ^/eJ^j(T*) models spontaneous emission, but for simplicity, we will consider only 
the noiseless case ej = 0. Also, we will assume the media are linear, nonabsorbing and 
nondispersive and the coupling is dissipative so that E^ = Ej and K is real. 

The parameter n in the first equation above is a symmetry-breaking coefficient measuring 
the extent of the feedback from the second laser to the first. The case of symmetric coupling 
(/i = l,Tc.i = Tc_2?Tf,i = Tf.2) has been analyzed in []13|. In this paper, we focus instead on 
the case of unidirectional coupling, i.e., /i = 0. When /i = 0, the first laser is unaffected 
by the second, and hence we may regard it as a driver. Assuming that this driving laser 
is pumped above its lasing threshold (pi > ai), it is easy to show that its amplitude and 
phase velocity settle down to constant values. Specifically, the long-term state of the drive 
is given by Ei = ^Jpi/ai — 1 exp (i($io + ^^i^*)), where $10 is an arbitrary constant, with 
constant gain Gi = ai. 

We assume these forms for Ei and Gi in the rest of the analysis. Thus, although we 
have formulated our study in terms of one solid-state laser driving another, the arguments 
and results we present hold for more general situations, such as a solid-state laser subjected 
to optical injection by other sources. 

In a typical experimental setup using Nd:YAG lasers |T^, Xf j and r^j are both positive 



and are of the orders 10"'* s and 10~^° s respectively. Thus we have two vastly different 
time scales in the system. The detuning ui — uj2 and the coupling K are control parameters. 
The detuning has values typically of order 10^ Hz, while the coupling can be varied over 
several orders of magnitude. In the symmetrically coupled system |T^, values of K ranging 



from 0(10 ^) to 0(10 ^) have been used. The pump pj and the loss aj are both positive, 
O(10~^), and their ratio Pj/otj is typically an 0(1) quantity. 



III. SCALING AND SYMMETRIES 

We scale the equations governing the response laser by introducing the following dimen- 
sionless quantities: 
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o = , B = , T = , F2 = — . 

a2Tf,2 "2 Tc,2 "2 

Here a is a stiffness parameter, typically O{10~^) in experiments, characterizing the vast 
difference in the time scales of the cavity round trip and fluorescence times in the response 
laser. The smallness of a will be important in the subsequent analysis. The parameter B 
is the dimensionless pump strength of the response laser; it often plays the role of a control 
parameter in what follows. The variables r and F2 represent dimensionless time and gain, 
respectively. 

Next we change variables by going into a reference frame rotating with the driver. Let 
X2 > be the amplitude of the complex field E2, defined by 

E2 = X2e'*^ 

and define $ by $ = $2 ^ '^'i; where $1 is the phase of the driving laser and $ ~ $2 " "^lo ~ 



{^i'^c,2/cx2)t for sufficiently large r. Assuming [T^l X2 7^ 0, we obtain the reduced system 

dX2 



dr 



(F2-l)X2 + fccos$ (3a) 

d^ ^ k . ^ , 
— = A sm $ (3b) 

dr X2 

^ = a{B-F2-F2Xi), (3c) 

where 

^ ^ (cu2 - iui)n,2 ^ ^ K^pilai - 1 
0:2 ' a2 ' 

Here, A is a dimensionless measure of the frequency detuning of the two lasers, and k can 
be interpreted as either a dimensionless coupling strength or injection amplitude. 

By choosing the phase difference $ as a variable, we have eliminated the explicit time 
dependence in the original system by rotating with the phase of the driver. In this rotating 
reference frame, a steady state now means a state in which the phase difference between 
the two lasers, and not the phase $2 of the second laser itself, is constant. Such a state is 
said to be phase-locked. The particular case in which $ = (which is possible if and only 
if A = 0) is called the coherent or in-phase state. 

The analysis of the reduced system (^ will occupy most of this article. As Eq. (|^) is 
invariant under 

A ^ -A, $ ^ (5) 

as well as under 

k -k, $ ^ TT + $, (6) 

we will assume from now on that fc, A > 0. Also, a, i? > by definition. 

There is a slight catch that one should be aware of. Although we can assume k > 
without any loss of mathematical generality, there can still be physical consequences. For 
example, (^ allows us to change the sign of k, but at the cost of transforming an in-phase 
solution to an antiphase one. In fact, for certain systems of coupled lasers, k can be negative 

m 
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IV. SPECIAL CASES 



The special case k = with general A is trivial: the driver and the response laser are 
decoupled. The reduced system (|) has a global attractor which is typically a periodic 
orbit, corresponding to phase drift between the laser and the drive. In degenerate cases, the 
attractor can be a fixed point for (||), e.g., if there is also no detuning, or if the pump for 
the response laser lies below the lasing threshold. 

A more important special case is that of zero detuning: A = with k ^ 0. Then $ = 
(in-phase) and $ = vr (antiphase) are invariant manifolds |]I8|, the former attracting all 



initial conditions except those on the latter, which is repelling. Hence, so far as long-time 
behaviors are concerned, we can confine our attention to the fiow restricted to the manifold 
$ = 0. 

On $ = 0, Eq. simplifies to the planar vector field 

"^^^ {F, - l)X, + k, 

a{B-F2-F^Xl). (7) 



dr 
dr 



The nuUclines IgU are 



1 - ^ (8a) 



(8b) 



1 + xi 



Since Eq. (|8a|) represents an increasing function of X2 while Eq. (|^) represents a decreasing 
one, there can be at most one real intersection. On the other hand, we may combine the 
two equations into a cubic polynomial equation in F2: 

F2[{F2 - If + A;2] - B{F2 - 1)' = 

and see that there is at least one real solution in the interval (0, min(l, 5)), by the in- 
termediate value theorem. Hence there is precisely one solution for F2 in (0, min(l, 5)), 
corresponding to a unique fixed point. An application of Dulac's criterion and the Poincare- 
Bendixson theorem to Eq. then reveals that this fixed point is actually globally stable 
with respect to perturbations on the X2F2 plane. Hence, in physical terms, a globally stable 
in-phase state exists if and only if there is no relative detuning between the two lasers. 

Having exhausted the possibilities for these special cases, we will assume A;, A 7^ from 
now on. 



V. LOCKED STATES 

A. Existence and stability for general parameter values 

A locked state of the laser system corresponds to a fixed point of Eq. (^, which must 
satisfy 
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(F2-l)X2 + fccos$ = 0, 
k 

A - — sin $ = 0, 
B-F2- F2XI = 0. (9) 

In the limiting case B = 0, (i.e., the response laser is not pumped), Eq. can be readily 
solved to yield 

k 

F2 = 0, X2= , $ = tan~^ A. 

The eigenvalues of the Jacobian at this fixed point are 

k^ 

Ai = -a(l + Y^^), A2,A3 = -l±zA, 

all having negative real parts. This means, as $ is the phase difference between the two 
lasers, that the response laser is passively driven into stable periodic motion 

E2 = -7=^== exp (z(tan-i A + $10 + ^^r)) 
V 1 + A^ 0-2 

with a constant phase-lead tan~^ A ahead of the driver, even though it is not pumped. 
For i? > so that F2 7^ at a fixed point, we may eliminate $ from Eq. (j^) and get 

Xl = ^- 1, (10a) 



Xl 



k^ 

(F2-l)2 + A2- 



(10b) 



Consider the graphs of X2 vs F2 corresponding to Eqs. ( |10a| , |10bD . The intermediate value 
theorem reveals that 

• if < i? < 1, there is only one fixed point, with F2 G (0, B) C (0, 1); 

• if 1 < B < Be, depending on parameter values, there may be one or three fixed points. 
There is always one in the interval F2 G (0, 1), and there may be a pair in F2 G (1, B); 

• if B > Be, again depending on parameter values, there may be one or three fixed 
points, all of which lie in F2 G (1, B); 

where 



Be = l + 



A2 



(11) 



represents a special value for the dimensionless pump strength. In Sec. |V B| we will show 
that Be, with a small correction, is the critical value at which locked states lose stability. 



7 



To analyze the local stabilities of these fixed points, we apply the Routh-Hurwitz criteria 
T9| , |20[| to the characteristic equation of Eq. and deduce that the fixed point is locally 



stable if and only if 



H2{F2) 



[B + 2)F| - 2F| - B{1 + A^) < 0, 



2a (F| - B 

2F2{F2 - 1 

2{F2 - 1)3 
2aB 



A1 < 0, 



F2 

+2(T^B 



-F2[{F2-l, 

-aB <0, 
2A'{F2 - 1) 

(F2-l)(F2-2)-2aF2(F2 

2a^B^ 



1) 



< 0, 



(12) 



where F2 is evaluated at the fixed point in question. When i/i, i^s < 0, the two inequalities 
H2 < and H2' < are equivalent and so we only need to check one of them. 

For B < Be, we can verify by direct substitution into (|I2|) that the fixed point with F2 < 1 
is stable. As we increase B through B^, this fixed point crosses -F2 = 1. The stabilities for 
the fixed points with F2 > 1 are not so clear. In fact, a.t B = Be we can solve the fixed-point 
equation ([T0|) exactly and find that either F2 = 1 or 



2A2 



± 



4A4 



A2 - P, 



where the last pair exist if and only if fc^/4A^ > A^ + A;^. Direct substitution into ([T^) reveals 
that the fixed point with F2 = 1 is stable (while the stabilities of the other two fixed points 
remain unclear) and hence there is no bifurcation at i? = Be- Hence, the condition B < B^ 
is sufficient but not necessary for a stable locked state to exist. In fact, for a = 0(1), we 
have found numerical examples of stable locked states even when B > Be- However, these 
examples are not of much physical interest, given that cr <^ 1 for real lasers. In the next 
subsection, we will show that for cr ^ 1, the locking condition B < Be is indeed tight, with 
a correction term that approaches zero as a does. 



B. Location and stability of locked states for fi ^ 1 

From now on we will assume cr -C 1. For definiteness, consider the scaling k,A = 0((t'^), 
and write 

k = Ka\ A = da". 

For a > 0, the graph of Eq. ( |10b| ) has a narrow peak, with 0((j") width, at F2 = 1, and 
we can show that there are three fixed points of Eq. (|^) if i? < B^, while there is only 
one fixed point if B > Be and B — Be ^ cx^". As we increase B, somewhere in the region 
Be < B < Be + 0{a'^"'), two of the fixed points collide. To better understand what is going 
on, we will employ the Routh-Hurwitz criteria ([T2|). 

For i? > 1, we can find the location of the fixed points perturbatively if a > 0. There is 
always a fixed point with 
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(13) 



For 1 < B < Be, there are a pair of fixed points witli 



\iB-l) 



2a\ 



(14) 



Direct substitution into (12) shows that only the fixed point with F2 = F* is stable. So 
we will concentrate on this fixed point and study the mechanism through which it loses 
stability. We expect that it will undergo a bifurcation in B^ < B < B^. + 0((T^"). 

Bifurcation occurs if one of the functions in (|12]) vanishes [^. Specifically, the Routh- 
Hurwitz condition Hi = corresponds to a zero eigenvalue, and is therefore generically 
associated with a saddle-node bifurcation, whereas = indicates pure imaginary eigen- 
values and a Hopf bifurcation. (These results follow immediately from the observation that 
the characteristic polynomial is cubic with real coefficients.) 

We solve (0) in conjunction with the fixed-point equation (p!0D by the method of domi- 
nant balance El. Let B = Bc + e. Then Hi{F2) = if 



e = El 



1 + 



4(5e - 1) 

a^-S^B, 



+ 0{a 



4:a\ 



2{B, - 1) 



+ 0(a 



4:a\ 



(15) 



and 1/3(^2) = if 



e 
F2 



2-2a 



1 + a^-^'^BeiB, 
for a < 1/2 



B^^{B,-l)/S' + 0{a 



3-4a\ 



(16) 



There is no self-consistent solution, with < e -C 1, for H^{F2) = if a > 1/2. If a < 
1/2, then <^ £1 and so becomes zero before Hi does, indicating a Hopf bifurcation. 
(Whether it is sub- or supercritical is undetermined at this stage in the calculation. But 
see Sec. |V11 B| .) If a > 1/2, then Hi becomes zero at e = Ei, indicating a zero-eigenvalue 
which, given the absence of any appropriate symmetry, is expected to be a 
[1^. We have confirmed these bifurcation scenarios numerically. 



bifurcation 

saddle-node bifurcation 
using the bifurcation package AUTO p4 



The argument above leaves a gap at a = 1/2, in which case Hi and H3 both vanish when 
e = 0{a), corresponding to a bifurcation of higher codimension, and potentially leading to 
a complicated outcome. We will look into this case separately in Sec. |VIII| . But as far as 
the fixed point is concerned, it loses stability at B = Bc + 0(cr). Hence, no matter whether 
the scaling exponent a is greater than, less than, or equal to 1/2, we can conclude that B^ 
is a tight estimate of the unlocking threshold, with an error at most of 0{a). 

To sum up, for cr <^ 1 we have shown that: 

• if B < Be + e, there is one locally stable fixed point. This corresponds to a stable 
injection-locked state, with a nonzero phase difference unless there is no detuning. 
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ii B > Bc + e, then there may be one or three fixed points, all in the region F2 G (1, B), 
and they are all saddles. Hence there is no stable locked state in this case. 



Here, < e = min(ei,e3) < 0{a), with 61,63 given by Eqs. (p^Sj) and (|1^), and B = Be + e 
is the parameter value at which the unlocking transition occurs. 

Taking advantage of the symmetries (P) and (^, we may rewrite the locking condition 
B < Be + e in terms of the physical parameters defined in Eq. (0) : 




P2 

- 1 > \UJI - UJ2\Tc,2\ (17) 

V "2 

with < e < 0(rc,2/tt2Tf,2)- In physical terms, stable locking occurs if the coupling K 
and the injection intensity (controlled by the pump strength pi) are sufficiently strong to 
overcome the relative detuning between the drive and the response laser. 

It is also noteworthy that the right-hand side of the locking condition depends on the 
pump strength p2 of the response laser. In applications involving injection locking, a laser 
with a low power output but very accurate frequency is often used to control a strong but 
sloppy one. The criterion ([TTD indicates that locking is more difficult if the response laser 
has a stronger pump p2. 

C. Discussion 

The condition for locking is commonly derived in the limit of small injection intensity 
(which equals pi/ai — 1 in our notation). This derivation relies on the slowly varying envelope 
approximation, along with the further assumptions that the intensity of the response laser 
is at the free-running output level and that the gain has saturated. Then the system of 
governing equations is reduced to a phase model equivalent to Adler's equation |]TT|. The 



traditional result for the locking condition is, in terms of our symbols, A < A;/vJ2 where I2 
is the free- running intensity of the driven laser. Since I2 = B — 1, this condition reduces to 
our condition B < Be- 

We have re-derived this well known result in a more general context. In particular, we 
have shown that a sufficient condition for locking is B < Be, even when k is not small and 
the phase model no longer holds. [Recall that k is the dimensionless product of the coupling 
strength and the injection amplitude; see Eq. (^.j Moreover, our approach allows us to 
calculate the correction terms in the formulas for the locking threshold, namely Eqs. (plSf ) 
and (|TB|), in the limit of small a = Tc^2/(tt2Tf,2)- 

To see why a simple phase model can give the correct locking condition, consider the 
least-stable eigendirection of the locked state. Suppose the coupling is weak, in the sense 
that k = 0(A) <^ a^l"^ . At the bifurcation point, the soft mode (the eigenvector associated 
with the zero eigenvalue) is given by 

(X2, F2) = (-a^ ^) + O(a-). 

For small a, the first and third components of this eigenvector are small, indicating that the 
soft mode lies nearly along the phase direction. Hence when the coupling is weak, the phase 
direction alone determines the locking condition, to a good approximation. 
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More generally, one can ask "When can the full model be replaced by a phase model?" . 
We will show in the following section that this reduction is valid when k and A are small. 



VI. WEAK COUPLING 

If B > Be + e, there is no stable fixed point of Eq. (|]), so the response laser cannot 
lock to the drive. The desynchronized dynamics then strongly depend on the relative sizes 
of the parameters k, A and a. For definiteness we will continue to assume that cr ^ 1 and 
k,A = Oia"). 

In this section, we focus on the weak coupling case a > 1/2, and show that unlocking 
occurs via a saddle-node infinite-period bifurcation. All solutions are then attracted to a 
stable periodic solution for the reduced system; this periodic solution corresponds physically 
to a phase- drifting state. 

Before we go on, a warning is in order: we will consider only lowest order effects in this 
section for simplicity. In particular, we will not pick up the e corrections [Eqs. ([T5|) and 
(|I^)] for the locking condition. 

Suppose that k = 0(A) <^ cr^^"^, or equivalently, 




\UJI - ^0-2 Tc,2 < 



rf,2 



in terms of the physical parameters in Eq. (Q). In this limit, we can reduce the full model 
Eq. @ to a phase model by the following argument. Let u = F2 — 1 and introduce the 
self-consistent assumption X2 = 0(1) (which we find to be justified by numerics). Then we 
observe that if m 3> A;, the term A;cos$ in Eq. (|3a|) is negligible and so the X2,u dynamics 
decouples from that of $ and constitutes a two-dimensional system. Next, by a phase plane 
analysis and diagonalizing the Jacobian we can show that X2 —>■ y/B — 1 and u decreases 
toward 0, with an 0(cr^/^) spiralling frequency and an 0{a) relaxation rate. So after an 
0{a~^) transient time, u becomes comparable to k and we can no longer neglect the term 
k cos $ in Eq. ( |5a] ) . 

As X2 = 0(1), Eq. (|3E| ) indicates that $ evolves with a characteristic rate of size 0{A,k) 
. If this characteristic rate is much smaller than that for the spiralling X2, u dynamics, which 
is true if A, A; ^ cr^^^i we can take fccos'l' to be a constant with an 0(A;^) error, and deduce 
that X2,'U relax to their equilibrium values with an 0(0"^^^) characteristic rate. 

Assuming this relaxation has happened, Eq. (^) becomes 

which is recognized as Adler's equation [Eq.(0)]. The important feature is that $ either 
approaches a constant or is strictly increasing, with $ mod 27r being periodic with period 
In I ^ — k"^/ (B — 1), depending on whether B is smaller than or larger than 1 + fc^/A^, 
with an 0{k'^) error. 

So, for large r, the unlocked dynamics are given by 
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kB cos^fr) 



2(5-1) ^ 

^2~l-^^^^ + 0(A;2), (19) 

where $(t) is governed by Eq. (|I8D. Hence X2(r) and F2(t) oscillate periodically, with 0{k) 
amplitudes and the same Fourier spectra as cos$., up to overall multiplicative constants. 
This last feature supplements the observation, noted in |^, p. 1153], that the laser intensity 
"acquires distortion at higher harmonics" as the locking threshold is approached. In fact, 
by insisting on a phase model, Siegman |]^, p. 1148] analyzed the dynamics of a laser with 
external injection outside the locking regime. Now we have justified this insistence by 
showing that in an appropriate regime in the parameter space, namely k = 0(A) <C cr^/^, the 
laser system indeed approaches a state with X2 — y/ B — 1, F2 ~ 1, starting from arbitrary 
initial conditions. 

The saddle- node bifurcation found in Sec. [VB| for a > 1/2 can now be identified more 



specifically as a saddle-node infinite-period bifurcation p. For B > Be + e, the phase 
increases strictly with a nonuniform speed, with a bottleneck near $ = sm^^{A\^B — 1/k). 
These behaviors are illustrated in Fig. |l|. 

The argument above clearly does not hold if A;, A = 0(o"^/^), as then the phase dynamics 
has a characteristic rate comparable to that of the spiralling motion on the UX2 cross- 
sections. We will address this case in Sec. IVlll. 



VII. STRONG COUPLING 

We turn now to the case of strong coupling: k = 0(A) ^ cr^/^. The main result is 
that unlocking occurs via a supercritical Hopf bifurcation, leading to a globally attracting 
limit cycle for Eq. (§) when B > Be + e. We will investigate in Sec. |VII A| the case k = 
0(A) = 0(1). Then in Sec. |V11 C we use a rescaling argument to show that this scenario 
also subsumes the ostensibly more general case k = 0(A) ^ o"^/^. 

A. Two-timing calculations 

Assume that a ^ 1, and k,A = 0(1). By defining 

E = Xae'*, 

which represents the complex electric field of the response laser as observed in a frame 
rotating with the driver, we can combine Eqs. (^aj) and (|3bD into 



^ = [(F2-1)+«A]E + A;. (20) 
ar 

For 0" -C 1, F2 is slowly varying and so may be taken to be a constant except over long time 
scales of order 0(l/cr). If F2 were really a constant, the exact solution to Eq. (|20D could be 
found as 
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E = Eeexp([(F2 - 1) +^A]r) - (21) 

where Ec is a complex constant. Hence X2 = \E\ would decay exponentially to 
k/ \J [F2 — lY + if F2 < 1 and grow exponentially to infinity if F2 > 1. In the first 
scenario where F2 < 1, we have 

ii B > Be and F2 — 1. So F2 would grow for F2 < 1. In the second case where F2 > 1, 
dF2/dT ~ — (TF2XI < and so F2 would decay. Hence we expect that, after initial transient 
instabilities for r ^ 0(l/cr), the system will settle down to its eventual fate, with F2 ~ 1 
and E ~ £'cexp(iAr) + ik/A for all sufficiently large r, and this is indeed what we find in 
our numerics (not shown). 

Now we investigate the structure of the attractor of Eq. (|^) by two-timing [^]. We first 
define a slow time T = err. Abusing notation, we have 

d df df 

where / denotes either F2 or E. Assuming transients have already decayed, we let F2 = 
1 + aF^^^ + 0{a^) and E = E^o) + aE^^) + 0{a^). Then we find 

^(0) ^ ^(0)(^y^)gi[9(0)(T)+Ar] ^ 



Fi'^ = f{T) + {B-l-[A^'\T)r~^jT 

2k 
A2 



+4^AW(r) cos(^^(°)(r) + Ar). 



To remove the secularity so that Fa^"* remains bounded (which is necessary for F2 to stay 
close to 1), we need 



= \Ib -{l + ^) = ^B- Be. (22) 



This is well defined if and only if 5 > i?c, which is, to lowest order, the condition for the fixed 
point of @j to be unstable (Sec. 0). Substituting the results for E^^^ and ^2^^ into the F^^-* 
equation and suppressing secularity so that it is possible to have E ~ Fcexp(2Ar) + ik/A, 
we get 

/ = 0, = ^T + = a^T + 

for some constant ?/'^*^\ Hence, for time r = 0(l/cr), and assuming the system is already on 
the attractor, we find 
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E 



ylWexp(z[(A + a^)r + ^W]) + | 



1 + 2o--J^v4W cos ((A + a^)T 



hO(a) 



(23) 



where = y/B - and is an arbitrary phase constant. These resuhs have been 
compared with numerics and show good agreement. Note that the result is singular if A ^ 0. 
A similar singular limit has been observed in a model of a CO2 laser p5|. 



B. Discussion 



1. Physical interpretation 
An inspection of E in Eq. (^) reveals that 

• if A^^^ < k/A, the phase difference $ oscillates between 7i/2 ± tan~^(Ay4'^*^Y^)- The 
response laser is said to be phase-trapped to the drive ||2^ : both lasers have the same 
average frequencies but their relative phase varies periodically. In other words, they 
are frequency-locked but not phase-locked. 

• if A^^^ > k/ A, then $ increases monotonically at an almost uniform rate, corresponding 
to a state of phase drift. 

These different behaviors, together with the regimes in parameter space in which they occur, 
are depicted in Fig. 0. This figure should be compared with Fig. |I], which shows that for 
k = 0(A) <^ 0"^/^, $ changes directly from locking to drifting at the critical pump value 
B = Bc + (o"^"), with no intervening possibility of phase trapping. 

Figure ^ illustrates the accuracy of our analytical approximations. For a given value of 
k, we compute the time series of the intensity / = X|. The local minima and maxima of 
the intensity are plotted in Fig. ^ as a function of k. For k > 2, the intensity is constant 
because the system has a stable fixed point corresponding to a locked state; hence only a 
single-valued branch of data is seen. As k decreases through the unlocking threshold, a limit 
cycle is born, causing two branches to bifurcate continuously from of the locked state, as 
shown in Fig. ^. This splitting is what one would expect for a supercritical Hopf bifurcation 
- the intensity now oscillates sinusoidally, so the lower branch corresponds to local minima of 
the intensity time series, and the upper branch corresponds to maxima. The curves passing 
through the data are analytical predictions given by Eq. (pSj). 

While the solution ( p3|) is periodic in the reference frame rotating with the driver, back 
in the laboratory frame the solution involves two frequencies. The phase of the receiver laser 
is then 

$2 = $1 + $ 

~ ^10 H -T + arg E 

0L2 

and is in general quasiper iodic. 
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2. Hopf bifurcation at B = Be + 0{a'^) 



The calculations in Sec. |V11 A] , and in particular Eq. (^), are correct only to the lowest 
order in a. If we proceed to higher orders, we find 



A 



From this and Eq. ( pSj) we see that the radius of the stable limit cycle scales as \/B — -Bbifn, 
where -Bbifn = -Be + 0{a'^) is the bifurcation value at which the limit cycle is born [compare 
with Eq. (|TBp], and the frequency A + ak'^/A^ is an 0(1) quantity. These results strongly 
hint that this is a supercritical Hopf bifurcation 0. Indeed, we can prove that this is the 
case for a -C 1 and k, A = 0(1), as follows. 

For B close to -Bbifn, we may take e = B — -Bbifn as a small parameter. By perturbing 
first with respect to a and then with respect to £, we find that the three eigenvalues at the 
relevant fixed point are 

Ai = ^(A + ^) + ^[^^ - - ^A)] + 0{a\ e% 



A2 — Ai, 

A3 = --5e + |^ + 0(a^.^), 

where an overbar denotes a complex conjugate. As e goes through zero, Ai and A2 cross the 
imaginary axis with nonzero speed, while A3 remains real and negative. Thus, by Hopf's 
theorem [^, a limit cycle is present on one side of the bifurcation point. Moreover, this 



limit cycle is centered at the fixed point, and has angular frequency given by the imaginary 
part of the conjugate pair of the eigenvalues. All these results are in agreement with Eq. 
(p3|), and we have also confirmed them numerically. 

To decide whether the bifurcation is sub- or supercritical, we use the Poincare-Lindstedt 
method to seek a periodic solution. We find that such a solution, be it stable or not, can 
exist only \iB > -Bbifn, i-e., when the fixed point is unstable, indicating that the bifurcation is 
supercritical. The periodic solution thus obtained is, as expected, the same as that given by 
Eq. (53). However, while the Poincare-Lindstedt method tells us nothing about the stability 



of the solution, it has the merit that it is uniformly valid for all time, thus confirming the 
speculation, based on numerics and the two-timing analysis in Sec. |VII A| , that there is 
indeed a periodic solution. 



C. Rescaling for cj < 1 with k = 0(A) > o^l'^ 

Having considered the case /c, A = 0(1), we will now demonstrate that as long as /c = 
0(A) ^ 0((T^/^), the lowest-order dynamics remain the same as before. Thus, any coupling 
and detuning that satisfy these milder bounds are still sufficiently strong to preserve the 
qualitative features seen earlier. In this sense, the condition k = 0(A) ^ 0(cr^/^) establishes 
the demarcation line for the regime of "strong" coupling and detuning. 
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By defining u via F2 = 1 + au, we may rewrite Eq. (j^) as 

dX2 



dr 



auX2 + k cos $, (24a) 



— = A-— sin<l>, (24b) 
dr X2 

'tL = B-{l + au){l + Xl). (24c) 

For A,k = 0(1), we observe in numerics that F2 = 1 + 0{a) so that u = 0(1). More 
generally, if fc, A = 0((j"), we may write A = Sa°',k = na"', with S,k = 0(1), and rescale 
time by setting t = a'^T. Then Eq. (p^D becomes 



dX-i l-2a 

dt 



ct^-^"mX2 + Kcos*, (25a) 



f = .-^s,„t. (25b) 



2 



^ = 5_(l + ^i-^)(l+X|). (25c) 
at 

For a < 1/2, our numerical simulations indicate that -F2 = 1 + {a^-") so that u = 0(1). 
Hence we can invoke a self-consistent argument to identify, to lowest orders, Eq. (pS]) with 
Eq. (|2^) , with a^"^", t, 5, k replacing a, r. A, /c respectively, since the terms au in Eq. ( P4c|) 
and a^~°'u in Eq. ( |25c| ) do not affect the lowest order phenomena and are therefore negligible. 
As 6,K = 0(1), we can apply the same two-timing analysis as in Sec. [VII A], but with the 



time scales t and ^"t, and get results similar to those found earlier. In terms of the 
variables used before rescaling, the lowest order terms are found to be 

Xse** = E = exp (z[(A + a-^)r + V^^")]) 

A'^ 

+ | + 0(a^--), 

F2 = l + 2a A^(o) ((A + a^)r + ^W) 

+0((t'~'"'^-^"), (26) 



where, as before, = y/B - = ^B-l- k^ / A"^ and ^/'(o) is an arbitrary phase con- 
stant. This result is the same as Eq. (^) to lowest order, exhibiting periodic motion born 
out of a supercritical Hopf bifurcation as B increases through B^. from below. One slight 
difference is that the next higher-order terms will be 0((T^~^") instead of 0(cr). Also, corre- 
sponding behaviors in this new system will occur on a slower time scale, as t = a°'T <^ r, if 
a > 0. 
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Clearly the argument above breaks down if a = 1/2, i.e., ii k = 0(A) = 0(cr^/^), 
as then cr^~^" in Eq. (|25a{) is no longer small as assumed by the perturbative treatment. 



Another way of putting this is, the two time scales, which are 0(1) and 0{a~^~^'^"') in t, 
collapse into one as a — >■ 1/2. In the next section, we turn our attention to this distinguished 
regime k = 0(A) = 0(cr^/^). For symmetrically coupled lasers, this is the regime where 
subharmonic resonance, amplitude instabilities, and chaos were discovered [0,0 . 



VIII. INTERMEDIATE COUPLING 
A. Overview 



The cases studied earlier allow us to anticipate some aspects of the possible unlocking 
behavior for intermediate coupling. From Sees. 0, |V|, and |V11| , we know that for weak 



coupling (a > 1/2), the fixed point corresponding to the stable locked state disappears in 
a saddle- node infinite-period bifurcation when B — = 0{a'^°'), while for strong coupling 
(a < 1/2), it loses stability in a supercritical Hopf bifurcation when B — B^ = 0((T^~^"). 
Hence, by sandwiching, we expect that for intermediate coupling (a = 1/2), the stable fixed 
point will lose stability when B ~ Be = 0(cr), probably in a codimension-two bifurcation 
combining the features of a saddle-node infinite-period bifurcation and a supercritical Hopf 
bifurcation. 

For B slightly above the bifurcation value, we expect the dynamics to be a combination 
of the two cases studied previously. The long-time behavior should be a spiralling motion 
that slowly passes through the bottleneck near 



X2 = VB^, $ = sin-i(^^^-^^), F2 = l, 

together with a fast re-injection roughly along the circle 

X2 = VB^, F2 = 1. 

Meanwhile, there is a global relaxation towards this attractor. The interplay among these 
three mechanisms can lead to complicated behavior, as we will see below. 

In contrast, in the case of symmetrical coupling (/i = 1 in Eq. (^), the phase difference 
decouples from the intensity and the gain, if one also assumes that both lasers have equal 
gains and intensities at all times |T^. This decoupling of the phase dynamics leaves only 



two mechanisms: a global relaxation toward the attractor, and a phase advancement (re- 
injection) along the attractor. In this sense, the unidirectionally coupled case might be prone 
to greater dynamical complexity. 

B. Numerics 

In the intermediate coupling regime k = 0{A) = 0(0-1/2)^ becomes 
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^ , 

— ; — = UX2 + K COS 

at 

| = B-(l + a"=«)(l + X|), (27) 

with K ~ 5 = 0(1) ^ cr^/^. In rectangular coordinates, with x = X2Cos$,?/ = X2sin$, 
this system can be rewritten as 

dx ^ 

— = ux — oy + K, 
dt 

— = uy + dx, 

^ = B-{l + a'/\){l + x' + y'). (28) 

To probe the dynamics, we will vary k and hold the other parameters fixed. This may 
seem unnatural, given that we have been using 5 as a control parameter so far. However, 
it turns out that if we vary k, instead of B, the features we wish to emphasize will stand out 
more clearly. In experiments, both parameters are easily tunable. 

Figure | illustrates the complicated dynamics that occur in this system. The plot is an 
orbit diagram, shown in the same format as Fig. ^ The local minima and maxima of the 
intensity's time series are shown as the coupling k increases toward the locking threshold. 
For the particular parameters chosen in the simulation, the locking threshold is at Kc = 
0.9986 = 6^/B — 1 + 0{a). Of course, for k > Kc, the results are simple: there is a single 
branch of constant intensity, corresponding to the locked state. But for k < we see an 
intricate pattern of local maxima and minima. 

The most striking feature occurs for 0.83 < k < Kc, where the diagram consists of a 
sequence of similar patterns that are scaled down as n increases. Figure ^(b) shows the 
region close to k = Kc in greater detail. The self-similar structure appears to persist all 
the way to the unlocking threshold at k = Hc- This scenario is not the standard period- 
doubling route to chaos, nor any of the other familiar bifurcation cascades. We do not fully 
understand what is happening here. Some first steps toward understanding this remarkable 



structure will be presented in Sec. IX. 



We are mainly interested in the self-similar structure near the unlocking threshold, but 
some other features of Fig. ^ deserve comment. An unsuspicious look at Fig. ^(a) might 
suggest that the attractor is quasiperiodic or chaotic for 0.52 < k < 0.58, as we see a 
complex set of data there. But this conclusion is premature, as a highly-looped periodic 
attractor with many turning points could also generate such a picture, thanks to the many 
local minima and maxima in its corresponding time series. 

Figure ^ shows that, in fact, this is exactly what is happening for many values of k, in 
the range 0.52 < k < 0.58. The xy projections of the attractors in phase space are plotted 
using DSTOOL 1 2^. The orbit in Fig. |^(a) appears to be quasiperiodic, while those in Figs. 



^(b)-(e) appear periodic. These figures also suggest that several kinds of of bifurcations are 
taking place. The trio Figs. ||(c)-(e) indicates a period-doubling sequence, although we do 
not see the full period-doubling cascade to chaos. 
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Hence, to distinguish whether the long-term behavior is periodic or not for a given set of 
parameters, a more careful treatment, such as plotting the power spectrum, is needed. We 
have computed the largest Lyapunov exponent as a function of k, with other parameters 
fixed, and find that it is always zero, up to the accuracy of the numerical algorithms. This 
result suggests that there is probably no chaos in the intermediate coupling regime considered 
here. Perhaps a more prudent way to put it is: chaos is not widespread for intermediate 
coupling, if it occurs at all. 

There are many other interesting attractors and bifurcations that could be discussed 
here, but we prefer to skip them so that we can focus our attention on the self-similar 
picture near the unlocking threshold. 



IX. A SINGULAR LIMIT 

We now consider a simpler system that exhibits the same cascade found in Fig. ^. Let 
cr = in the reduced system (|27|) , so that it becomes 

— ; — = UX2 + K COS <P, 

at 

§ = B - (1 + A1). (29) 

Figure ^ shows that this new system still exhibits the striking pattern found earlier 
(although the precise positions of the bifurcations are not the same, and there are structures 
that come up in one case but not the other.) This robustness suggests that the self-similar 
cascade is not due to perturbative effects of a but is more generic. To check that the cascade 
is not just an artifact of the choice of variables plotted, we have also tried plotting local 
minima and maxima of x, y, u, and I = + against the parameter k. In all cases, we 
find similar patterns. We have also checked that the patterns keep recurring even closer to 
the unlocking threshold, specifically for k G [0.998, 1] and k G [0.9998, 1]. So we believe that 
the self-similarity is genuine. 

On the other hand, the limit cr = requires some caution on our part. The blown- 
up system ( |27| ) is obtained from Eq. (|^) by the rescaling t = a^^'^r, a^^'^u = F2 — 1, A = 



cr which becomes singular as cr 0. Hence Eq. ( p9D is not directly rele- 

vant physically, but it has the virtue that it leads to cleaner numerics, and brings out the 
bifurcation mechanism more clearly. 

One should also be aware that the singular limit displays some nongeneric dynamics 
because of its higher degree of symmetry than the original system. For instance, Eq. ([29| ) is 
invariant under the transformation 

T : t ^ —t, X —X, y ^ y, u ^ —u, (30) 



a "reversibility" symmetry that is not preserved by Eq. (|27| ) for cr 7^ 0. A corollary of 
this invariance is that if Q is an attractor, then the set Q = {{x,y,u)\{—x,y,—u) G Q} 
is a repellor. Results like this can be used for checking numerical convergence in tracing 
bifurcation diagrams. 
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A. Fixed-point analysis 



Since the transformation from Eq. (^) to Eq. (^) is singular, the fixed point analysis 
(both existence and stability) in Sec. |V|is not directly applicable to Eq. (pQ]). The necessary 
modification is straightforward, and the result is as follows. Fixed points exist if and only 
iil<B<Bc = l + They are given by p± = {±x*,y*, ±u*), where 



X 



u 




B-l){B,-B) 



B,-B 
B-1 ' 



Hence both fixed points recede to infinity as B 
B B-. 

As for stability, the Jacobian of Eq. ( p9|) is 



J 



u —5 X 
6 u y 
-2x -2y 



and they coalesce at (0,k/(5, 0) as 



At the fixed points, the characteristic equation simplifies to 



= A(A) = - 2uX^ + 



B - 1 



A + 2(5- l)(A-n), 



(31) 



where u = at the fixed point p±. Since /-«(— A) = — /u(A) and u*^_ = — u1, we can 
conclude that if A is an eigenvalue of p+, then —A is an eigenvalue of Hence we cannot 
have a saddle-node bifurcation ai B = Be- 

The eigenvalue equation (|3l|) , as it stands, is hard to solve. But some information can 
be obtained from the intermediate value theorem. Since /„(±oo) = ±oo, and /«-^(0) = 

±2(5 — l)5-y/ {Be — B)/{B — 1) is positive (negative) at p+ (p-), the intermediate value 
theorem implies that the fixed point p_ always has a positive eigenvalue and therefore 
cannot be stable. Whether it is a saddle or a repellor is not determined by this argument. 

To investigate the type of bifurcation at B = Be, we will consider Eq. (^) with B = 
Be — e, where e > is a small parameter, and k, and 6 are assumed to be 0(1) positive 
quantities. Then we have a regular perturbation problem and can find the eigenvalues to be 



Ai,2 = ±i^S^ + 2kV52 _ 

^3 = ^i:p^V^ + o(.) 



k6^ + 2k^ 



(32) 



The 



at p+. Hence p+ is an attractor. More specifically, it is an attracting "spiral-node", 
eigenvalues at p_ can be obtained by taking the negatives of the eigenvalues of p+. 

These results provide some local information about the bifurcation at i? = B^. As 
B B~, i.e., e —>■ 0~^, the two fixed points approach each other and collide at i? = 
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Be- Meanwhile, their attractiveness and repulsiveness get weaker and weaker, since all the 
eigenvalues approach the imaginary axis. The fixed point ai B = is linearly neutrally 
stable. 

As Eq. ([3^ ) indicates, locally the bifurcation is of codimension two, with a simple zero 
eigenvalue and a conjugate pair of pure imaginary eigenvalues. As noted in |jT8|, such 



bifurcations can be complex because of their global structures. In Eq. (^), as in Eq. (P4D, 
the re-injection provides such a global mechanism that leads to complicated dynamics. 

Since the stable fixed point is annihilated in a collision with an unstable object as -B — > 
i?c , we do not expect stable objects for B > B^to be in small neighborhoods of (0, 5{Bc — 
0), the position of the fixed point at i? = B^.. However, some form of intermittency 
might be expected, and as we show in the next section, this is indeed the case. A trajectory 
on the periodic attractors spends a long time in helical motion near the ghost of the fixed 
point. The same is true when cr > 0, indicating that this helical motion in the bottleneck is 
not caused by the degeneracy. 

Before we move on to numerics, it should be noted that this degenerate "attractor-repellor 
bifurcation" is a consequence of setting a to zero. For a > 0, a generic saddle-node bifurca- 
tion occurs at i? = B^. Also, for o" > 0, generically we do not have a genuine codimension-two 
bifurcation. But since both Hi and in (|T2p vanish at B^. + 0((t), meaning that all three 
eigenvalues are close to the imaginary axis when bifurcation occurs, the nonsingular system 
(^) may pick up some remnants of the bifurcation scenario of the singular system (|29|) , so 
that their orbit diagrams (Figs. ^ and exhibit similar features. This suggests that what- 
ever the bifurcation scenario is, it is stable with respect to perturbations in parameters, at 
least in an operational sense. 



B. Numerical studies of the self-similar bifurcation sequence 

In our numerical studies of the self-similar bifurcation sequence, we concentrate on the 
periodic windows. As before, we will vary k and keep other parameters fixed. Then for 
K = + e, with < e <^ 1 and = 5\/ B — 1, the fixed points are at p± = {±x*,y*, ±u*), 
where 

He 

u* = v^(l + 7^) + 0{e'/^). 



The eigenvalues at p+ are 



2 + 
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1. Helical structure of the periodic attractors 



Using AUTO, we have examined the phase-space geometry of the attractors in individual 
periodic windows. For definiteness, consider the yu projections of the attractors. Figure |^ 
shows a typical collection of snapshots as we sweep k, across a window. What happens in 
this window also happens in other windows, with some modifications in details. Specifically, 
we see n little helical loops per round trip at the left end of the nth window. (In this 
example, n = 7.) Most of these loops are located near (0, 1,0), where the stable fixed point 
is annihilated at k = 1. 

Intuitively, the trajectory is slowly funneled through the bottleneck caused by the ghost 
of the former stable fixed point. During this slow passage, the trajectory also spirals around 
at an 0(1) frequency, given by the imaginary part of the complex eigenvalues in Eq. (P3[). 
The combination of slow passage and 0(1) spiralling gives rise to the helical appearance of 
the trajectories. 

As we increase k in a periodic window, (n — 2) of these loops stay at roughly the same 
distance from the u axis, while the remaining two loops are stretched in the re-injection part 
of the flow. Moreover, the left side of the projection moves, with a stronger shear at the 
bottom, towards the u axis, eventually crossing the loops. Meanwhile, there is a clockwise 
rotation to bring about the structure we see as we reach the right-hand end of the window. 
The same trend occurs in all the periodic windows we have sampled, with 3 < n < 12 (we 
do not see periodic windows corresponding to n = 1, 2). If this "n loops in the nth window" 
scenario is correct, and if there really exist infinitely many such windows accumulating at the 
limiting value k = 1, we will have an attractor with infinitely many tight loops at k = 1~, 
just before its annihilation in an attractor-repellor bifurcation. This will then be reminiscent 
of the chaotic orbit born from the infinite sequence of period doubling in the logistic map. 

In a similar problem concerning a CO2 laser, Zimmermann et al. worked backwards 
and constructed a model by assembling a local spiral motion and a global re-injection. 
Numerically, for a fixed set of parameters, the flow is similar to what we see here. However, 
how this mathematical model is related to the original laser system is unclear. 



2. How periodicity ends 

In each periodic window, the unstable branch always has its nontrivial Floquet multipliers 
outside the unit circle, indicating a repellor instead of a saddle. When this branch collides 
with the stable branch, an attractor-repellor bifurcation of cycles occurs, bringing an end to 
the periodic window. In virtue of the invariance of Eq. (^) under the transformation T in 
(pop, it is no wonder that at the bifurcations, when the attractor and the repellor coalesce, 
this neutrally stable object is invariant under T' : a; — > — x, y —>■ y,u ^ —u. This explains 
why the projections onto the yu plane are always symmetrical with respect to u ^ —u [as 
seen in Figs, ^(a) and (d)], while the projections onto the xy plane are always symmetrical 
with respect to a; — > —x at the ends of the windows (not shown). 

As mentioned earlier, the attractor-repellor bifurcations are a consequence of the re- 
versibility symmetry of the cr = system. Generically, the periodic attractors are annihi- 
lated in saddle-node bifurcations and there is no T-symmetry, a strong indication that the 
mechanism underlying the cascade of bifurcations has nothing to do with this symmetry. 
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Recalling the results in Sec. [IX A| , we notice that while an attractor-repellor bifurcation 
of fixed points occurs as the locking threshold is approached from above {k — >■ l"*"), there 
is a simultaneous attractor-repellor bifurcation of cycles (presumably with infinitely many 
loops) as the threshold is approached from below (/t — > 1^). This double-sided aspect of 
the unlocking bifurcation also seems to be preserved for a > 0, suggesting that it too has 
nothing to do with symmetry. 



3. Trends in the periodic windows 

The periodic windows in Fig. ^ have been verified by auto to be isolas, i.e., they ter- 
minate at both ends as the periodic attractors are annihilated in collisions with unstable 
objects (probably unstable periodic orbits). However, we do not fully understand the mech- 
anism that causes such annihilations, although it seems to be some form of resonance of the 
spiralling motion and the global re-injection. Indeed, we have found that as we increase k in 
a periodic window, periodicity ends precisely when a trajectory executes an integer number 
of spiral loops upon one re-injection, i.e., if 

27r 27r 

— = n—, 34 

uji Us 

where n is a positive integer, and ui and us are the frequencies of the re-injection and the 
spiralling motion, respectively. To see this, we assume that the time needed for re-injection 
is dominated by the slow passage through the bottleneck. Then, using Eq. (|33|), we may 
rewrite Eq. (|^) as 



\3/2 



(54 + 2Ke2^ n 

to lowest order, where En = ^ — K,R,n, and KB.,n is where the nth periodic window ends at the 
right. Solving for e„, we get 

_ (^^ + 1 

S^{2K,r ^^^^ 

Moreover, at these parameter values, corresponding to the right ends of the periodic windows, 
the periods of the attractors are expected to be 

271 5^ + 2kc^ 1 27r(5 



Figures ^(a) and (b) compare these predicted scaling laws against numerics. For n > 5, 
Eqs. ([35| ) and (|36D agree with numerics to within 2%. 

We do not know what mechanism kills the periodic window as we decrease k. But drawing 
analogy to the mechanism at the other end, we expect that the positions of the left ends of 
the periodic windows scale as 

1 - /«L,n = 0(l/n=^), 
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so that the widths of the periodic windows scale as 

Wn = KR,n- HL,n = 0{l/n^), (37) 

and the periods of the attractors at the left ends of the periodic windows scale as 

TL,n = 0{n). (38) 

These scaling behaviors are verified in Figs. |Kc)-(d). 

Meanwhile, the sizes of the attractors, as measured by the L^-norms, remain 0(1) as 
K — s> 1~. This is expected as the L^-norm is dominated by the global re-injection, which 
persists in all periodic windows. 



X. OPEN QUESTIONS 

From a theoretical perspective, the most interesting open question concerns the mathe- 
matical mechanism underlying the self-similar cascade of bifurcations observed in both Eq. 
(p7|) and the simpler system (pQ]). Whatever the mechanism is, the heuristic arguments 



and numerical evidence presented in Sec. |I^ suggest that it must combine the features of a 
saddle-node infinite-period bifurcation and a supercritical Hopf bifurcation. As such, it may 
well arise in other scientific settings. Maybe it can even be detected experimentally. 

Although self-similarity itself is common in dynamical systems, and has been explained 
by renormalization group arguments in such contexts as period doubling, intermittency. 



and quasiperiodic breakdown |18[, it seems the cascade we see here falls into none of these 
categories. Rather, it is characterized by an infinite series of saddle-node bifurcations of 
cycles, accumulating at a finite parameter value corresponding to the locking threshold. 

Aside from these bifurcation issues, the dynamics of Eq. ( [25| ) is also interesting for a fixed 
set of parameters. Zimmermann et al. suggested that the periodic orbits with many 
small loops can be understood as the effects of a flow with helical local dynamics together 
with a global re-injection. The challenge now is to find a way to reduce the laser equations 
to a form where this conjectured phase space geometry becomes transparent. 

There are many other interesting avenues for future research. The dynamics of the 
unidirectionally forced system could be explored over a much broader range of parameter 
values, with k and A not necessarily of the same order. We have also neglected the effects of 



noise, a topic of great importance in technological applications of injection locking [|ri| ,p0|-|32| 



Another promising direction would be to study arrays of coupled lasers PB|-P^ driven by 
external injection, particularly in regimes where amplitude effects are important and the 
phase model approximation is not valid. Finally, we have restricted attention to drive 
signals of constant intensity and frequency, but for applications to optical communications 
40| , |4l[| , one needs to study how lasers respond to modulated drive signals, especially those 



carrying messages within them. 
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FIGURES 



FIG. 1. The dependence of the behaviors of the phase difference ^ on the pump to loss ratio B, 
for A; ~ A = 0((t"), a > 1/2. A, k are as defined in Eq. @. E = ^26** is the complex electric field. 
The stable (unstable) fixed point is represented by a "•" ("o"). $ is constant ii B < Be = l+fc^/A^ 
and increases strictly if B > Be, with a saddle- node infinite-period bifurcation at B = Be- The 
bifurcation value for B has 0(ct^") errors [Eq. (|l5|)] and the value of X2 has 0{a'^) errors [Eq. (|l9|)]. 
Contrast this with Fig. ^. 

FIG. 2. The dependence of the behaviors of the phase difference $ on B, for 
A; ~ A = 0{a"'),a < 1/2. Parameters, variables, and symbols are as defined in Fig. |^. $ is 
constant if B < B^, oscillates between 7r/2 ± tan"^(AA(°V^) if -Be < -B < 1 + 2k'^/A'^, and in- 
creases strictly if B > 1 + 2k'^ / A? . ^(o) is as defined in ([2^). A supercritical Hopf bifurcation 
occurs at = B^. The bifurcation value for B has 0(cr^~^") errors [Eq. (16)] and the value of X2 



has 0((T^ ^") errors [Eq. (pB])]. Contrast this with Fig. |l|, where phase trapping is impossible. 

FIG. 3. Orbit diagram, plotting the local minima and maxima of the intensity / = Xl vs 
/c, of Eq. (|3|), with a = 0.010, A = 2.0, /c G [0, 2.2],i? = 2.0. Dots represent numerical results and 
curves theoretical expectations [Eq. (p3|)]. In this and subsequent orbit diagrams (Figs. ^ and 
^), the numerical results are obtained by following the equation of motion with the subroutine 
STIFF in [W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes 
in Fortran: The Art of Scientific Computing (Cambridge University Press, New York, 1992), 
2nd ed.], which implements a fourth order Rosenbrock method, and then estimating the extremal 
values of the intensity with a second order interpolation. The results are corroborated by another 
subroutine, stifbs [ihid], which implements a semi-implicit extrapolation algorithm analogous to 
the Bulirsch-Stoer method. 



FIG. 4. Numerically computed orbit diagram of Eq. (^, with a = 0.0025, 5 = 1.0, B = 2.0. 
(a) K G [0.40,1.1]; (b) a blow-up in the range k G [0.98,1.0]. Contrast the complicated patterns 
with the simple structure in Fig. ^. Note in particular the self-similar bifurcation sequence that 
piles up at K = 0.9986 = + 0{a). 



FIG. 5. The xy projections of the attractors of Eq. (g^), with a = 0.0025,5 = 1.0, S = 2.0. 

(a) K = 0.78, a quasiperiodic attractor. We have followed the flow only for a time span of 1000 
units in order to reveal the structure, (b) k = 0.77, a periodic attractor with many local minima 
and maxima. This attractor remains a closed curve without blurring in a time span of 2000 units, 
strongly hinting its periodicity, (c) k = 0.70, a period-1 orbit, (d) k = 0.68, a period-2 orbit, (e) 
K = 0.67, a period-4 orbit. 

FIG. 6. Numerically computed orbit diagram of Eq. (29), with 5 = 1.0, B = 2.0. (a) k G [0, 1.1]; 

(b) K G [0.98, 1.0]. Compare with Fig. |. 
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FIG. 7. The yu projections of the attractors of Eq. (29) as we sweep across the n = 7 periodic 
window by increasing k, with S = 1.0, B = 2.0. (a) k = 0.91867; (b) k = 0.92500; (c) k = 0.92980; 
(d) K = 0.93200. The snapshots are not taken at evenly spaced values of k. Rather, representatives 
are chosen to show the deformations more clearly. 



FIG. 8. Log-log plots showing the trends among the periodic windows of Eq. ( [29| ) with 
5 = 1.0, B = 2.0, in the (a) positions of the right ends of the periodic windows, with the theoretical 
prediction ln(l — KR,_n) = ln(27/8) — 2 Inn; (b) periods of the attractors at the right ends of the 
periodic windows, with the theoretical prediction InTR^n = ln(27r/\/3) + Inn; (c) widths, with the 
fit In Wn = 1.5 — 3 Inn; and (d) periods of the attractors at the left ends of the periodic windows, 
with the fit InTL^n = 1-2 + Inn. Dots represent raw data with 3 < n < 12, solid lines represent 
theoretical predictions [Eqs. (]35| ) and (|3^)] with no fitting parameter, and dashed lines represent 
semi-theoretical predictions, where the slopes are chosen according to the theoretical scaling [Eqs. 
(37) and (|3^)], and the intercepts are chosen to fit the data. 
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